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1. Introduction 

The primordial state of matter called quark-gluon plasma (QGP) is expected to be realized in 
extremely hot and dense mediums, and a lot of experimental efforts have been made to produce 
such a state in heavy-ion collision experiments [|I]]. To understand QGP, theoretical studies by 
the first principle calculations of QCD at high temperature and density are important. At present, 
the lattice QCD simulation is the only systematic method to do so. Many important properties 
of finite temperature QCD have been studied by the Monte-Carlo simulations [Q]. The studies at 
finite density had been known to be difficult until recently. However, recent technical developments 
allow us to extract information on those in the low density region. In this report, we would like 
to review recent progress in lattice QCD at finite density. One of the most interesting studies is 
to investigate the phase structure of QCD at non-zero temperature and density. The QCD phase 
transition has been found to be crossover at zero density by simulations with staggered type quark 
actions [||]. We expect that the nature of the phase transition will change to be of first order in the 
high density region, and it is very important to find the critical point terminating a first order phase 
transition line, since the critical point is one of the most interesting features that may be discovered 
in heavy-ion collision experiments. The study of the equation of state (EoS) is also important. The 
numerical studies by lattice QCD simulations will be able to provide basic input for hydrodynamic 
calculations of the expansion of hot and dense matter generated in heavy-ion collisions. 

However, lattice QCD at non-zero density is known to have a serious problem. In a Monte 
Carlo simulation, we generate configurations of link variables {U^(x)} with the probability in 
proportion to the weight factor (detM) Nf e~ Sg and the state density of {U^ (x)}. Here, M is the quark 
matrix and S g is the gauge action. The expectation value of an operator ^[U^] is then evaluated by 
taking an average of over the generated configurations {U^(x)}. 

^Vconf. {U)i{x)} 

The quark matrix at zero density have the 75 Hermiticity and the Hermiticity guarantees that the 
quark determinant is real. However, the relation of the 75 Hermiticity changes to 

M\n q ) = Y5M(-^)y 5 . (1.2) 

at finite quark chemical potential (p. q ). Then, the quark determinant becomes complex except for 
\i q = 0; (detM(/i 9 ))* = detM(— \x q ) / detM(jU^). Because the Boltzmann weight must be real 
and positive in the Monte-Carlo method, we cannot perform a simulation at finite density directly 
A popular method to deal with QCD at finite \L q is the reweighting method. However, we will 
encounter another problem called the "sign problem" in the calculation at large \i q . The key point 
in the study of finite density QCD is to avoid this problem. 

A lot of progresses have been obtained in this field. The equation of state in the low density 
region was studied in |Q|5[ |(| [7|]. The sign problem is still one of the most important issues in the 
study of finite density lattice QCD. The nature of the sign problem was discussed in the random 
matrix model and the chiral perturbation theory []|, [J |Kj. Some trials to find the critical point 
at finite density were examined [11, 12, 13, [14]]. Moreover, a new algorithm based on stochastic 



quantization was proposed in [15]. The phase structure in the high density region was studied in 
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the strong coupling limit [ |16| [17]]. The hadronic fluctuations in the high temperature phase was also 
studied using an effective theory U18Q. Moreover, the equation of state by chiral fermion actions 
was discussed in the high temperature limit [ |l9[ ^Oft . The phase structure of two-color QCD, which 



is free from the sign problem, has been studied. (See e.g. [ |21| , [22[ , |23| ] for a review.) Among these 
topics, we want to focus on the equation of state and the critical point in the (T,/j, q ) plane in this 
report. We discuss the equation of state and hadronic fluctuations in Sec. ||[ Some attempts to find 
the critical point are discussed in Sec. |[ A summary is given in Sec. ||. 



2. Equation of state at finite density 

In order to extract unambiguous signals for the QCD phase transition from the heavy-ion col- 
lisions, quantitative calculations from the first principles of QCD are indispensable. In particular, 
studies of the equation of state (EoS) can provide basic input for the analysis of the experimental 
data. Many studies have been done at finite temperature (T) and zero chemical potentials (jj, q ) [0]. 
Also, recent developments of computational techniques enabled us to extend the study to small \i q . 

Several years ago, systematic simulations for the study of the EoS at finite density have been 
performed by the Bielefeld-Swansea Collaboration using p4-imploved staggered quark action with 



rather heavy quark masses Q24| , |25| , [26; ] . They found that the Taylor expansion method is useful for 
the EoS study in the low density region which is important for heavy-ion collisions. Moreover, they 
found large fluctuations in the quark number density at finite density. The temperature dependence 
of the quark number susceptibility % q , which corresponds to the fluctuation of the quark number, 
changes qualitatively when jx q becomes non-zero. For \i q = the susceptibility Xq/T 2 changes 
rapidly at the transition temperature but continues to increase monotonically. However, for \i q ^ 
the quark number susceptibility develops a pronounced peak at the transition temperature. Such a 
behavior suggests the existence of a critical point in the {T,ji q ) phase diagram. 

In this year, remarkable results were obtained by simulations near the physical quark mass 
point with improved staggered quark actions [Q, ||| q|. The MILC Collaboration and the RBC- 
Bielefeld Collaboration studied the isentropic equation of state, i.e. the EoS along trajectories of 
constant entropy per baryon number. There are also progresses in the study of fluctuations at fi- 
nite density. The RBC-Bielefeld Collaboration found that the enhancement of the quark number 
susceptibility becomes larger as the quark mass decreased Moreover, the WHOT-QCD Collab- 
oration performed simulations with a Wilson type quark action and studied the EoS at finite density 
P]. They calculated the quark number susceptibility and confirmed the large fluctuation at \i q ^ 0. 

2.1 Taylor expansion method 

The main problem in the study of QCD at finite density is that the Boltzmann weight is com- 
plex for \i q ^ 0. Because the Boltzmann weight must be real and positive if we want to generate 
configurations with the weight, the conventional Monte-Carlo method is not applicable at \i q ^ 0. 
One of the possible approaches to study the finite density QCD is performing a Taylor expansion 
of physical quantities in terms of \i q around \i q = and calculating the expansion coefficients by 
numerical simulations at \i q = [Q, 25, 26, 27, 28]. Because the simulations at \i q = is free from 



the complex weight problem, the expansion coefficients, i.e. derivatives of physical quantities with 



3 



Recent progress in lattice QCD at finite density 



Shinji Ejiri 



respect to H q /T, can be evaluated by a conventional Monte-Carlo simulation. The pressure (p) is 
obtained from the partition function (5°), 

4- = ^ln3f = a, (2.1) 

and the calculations of the derivatives the partition function are basic measurements in the QCD 
thermodynamics, since most of thermodynamic quantities are given by the derivatives of £2. 
We define the Taylor expansion coefficients as 

p__ y u,d,s (T) (Vu VY^y/M* uAs = d i+ j +k a 

T 4 i,tJ 1 ' 3 ' ^TJ \ t) \t) ' i\jlk\ dOiu/TydQia/Tydfa/T)' „ 

(2.2) 

Here, /v</,s are the chemical potentials for the u,d,s quarks, and Cq'qq(T) is the pressure at \i u = 
l±d = Ms = 0. The coefficient c?'j£(T) are computed by performing a simulation at \i q = 0. The 
explicit forms of the Taylor expansion coefficients are given in [Q, ^6| ]. We expect that QCD in the 
high temperature limit is described as free gas of quark and gluon and the \i q -dependence of p/T 4 
is given only through terms of }x q and \L q for the free gas. Therefore, the Taylor expansion may 
converge well in the high temperature region. 

For the calculation of pressure at \i q = 0, the integral method is commonly used. Using the 



thermodynamic relation Eq. (|2.1|), the pressure is computed as 



>-vk i **-3fi = -vk i *\-w)- <23) 

Here, 5i at is the lattice action and (• • •) is the thermal average with zero temperature contribution 
subtracted for the normalization of p. In multi-parameter cases such as full QCD, j8 should be 
generalized to the position vector in the coupling parameter space. The initial point of integration 
/3o is chosen in the low temperature phase from the condition p(fio) ~ 0. The derivatives of S\ at 
with respect to /3 and the quark mass are basically given by the Wilson loops and chiral condensate. 
The energy density is obtained from the following equation, 



3p 1 dhiZT 



T 4 VT 2 dT 




a dS ^\ _JP /dflat , (?4) 

, ( , i \ da / da \ dp 



where a is the lattice spacing, the lattice size is N 3 x N t , and j3 is the position vector in the coupling 
parameter space for full QCD, again. The density effect of (e — 3p)/T 4 can be estimated by a 
Taylor expansion. The coefficients are given by the derivatives of Eq. (0) with respect to jJ. q /T. 
The quark number density n u j s is calculated by 

n M ,,_ 1 dlniT _ d(p/T 4 ) 



and Eq. (2.2). We define the light quark number density as n q = n u + n c i. The susceptibilities of 
light Quark number (% q ) and strange quark number {% s ) are given by 

X q _ ( d + d \ n u + n d Xs_ = d(n s /T 3 ) ^ g) 



T 2 \d{pu/T) d{}Xd/T)) J 3 ' T 2 d(ix s /T) 
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Figure 1: Pressure (left) and energy density (right) vs. temperature along the lines of constant entropy per 
baryon number obtained by 2+1 flavor simulations with asqtad staggered fermion action [gj], 




Figure 2: (left) The ratio of pressure and energy density as a function of energy density on isentropic 
trajectories obtained with p4fat3 staggered fermion action for N t = 4 (filled) and 6 (open), (right) The ratio 
of %l and %2 vs. temperature for m K w 220MeV [|] and m, « 770MeV [|(|. 



These susceptibilities correspond to the fluctuations of the quark numbers. Moreover, the entropy 
density s is given by the thermodynamic relation, 

j_ _ e + p-Zf=u4^f n f n 7 x 

The chiral condensate is defined by the derivative of In with respect to the quark mass. 

2.2 Isentropic equation of state 

One of the most interesting results which have been obtained from heavy-ion collision exper- 
iments is that the experimental data is well-explained by a perfect fluid model without viscosity. 
This implies that a dense medium created in a heavy-ion collision expands without further gener- 
ation of entropy after thermalization. Therefore, it is important to calculate the EoS with keeping 
the entropy (S) per baryon number (Nb) constant for the analysis of the experimental data 

The MILC Collaboration and the RBC-Bielefeld Collaboration studied the isentropic equation 
of state by performing simulations near the physical quark mass point using improved staggered 
fermion actions. The isentropic expansion lines for matter created at RHIC, SPS and AGS ener- 
gies correspond to S/Nb ~ 300, S/Nb ~ 45, and S/Nb ~ 30, respectively. These values have been 
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obtained by comparing experimental results for yields of various hadron species with hadron abun- 
dances in a resonance gas [p0|]. Measuring the Taylor expansion coefficients of the pressure, energy 
density, baryon number and entropy by Monte-Carlo simulations, they found the isentropic trajec- 
tories in the (T ,pL q ,pL s ) parameter space with an additional constraint n s = 0, where \L U = \id = AV 
They then calculated the energy density and pressure along the isentropic trajectory using the Tay- 
lor expansion coefficients. 

The results of energy density and pressure computed by the MILC Collaboration are shown in 
Fig. [I] for each S/Nb @, 0]. The pion mass is about m n m 220MeV, which is close to the physical 
pion mass. They used the asqtad quark action and successfully reduced the discretization error in 
the EoS. The filled and open symbols in these figures are the results on lattices with a temporal 
extent N t = 6 and 4, respectively. The difference between them is found to be small. 

The RBC-Bielefeld Collaboration also studied the isentropic equation of state performing 2+1 
flavor simulations with m n 220MeV using the p4fat3 action The results of the pressure and 
energy density are consistent with the results by the MILC Collaboration. The right panel of Fig. || 
is the results of p/e plotted as a function of e. The open symbols are the results from N t = 6 and 
the filled symbols are from N t = 4. They found that the density dependence of p/e is small when 
they plot as a function of e. Also, the speed of sound in the dense medium, c 2 s = dp/de, can be 
calculated by measuring the slope of p/e in this figure. 

2.3 Radius of convergence and hadronic fluctuations 

Next, let us discuss the convergence radius of the Taylor series, p/T 4 = Lcf (/i B /r) ! , with 
cf = (l/3)(l/il)(Y,f= u d s d/d(iJ,f/T))'Q.\i luds= o. We expect that the crossover transition at low 
density changes to a first order phase transition at a critical value of /ig. If there is such a critical 
point, the Taylor series does not converge at the critical point. The simplest way to estimate the 
radius of convergence (p) is to calculate the ratio of the expansion coefficients. We define 

p = lim p n , p n = J\c B /c B 2 \. (2.8) 

7)^00 V r 

In the region of Hb/T < p, p/T 4 is finite. For the case of free quark gas expected in the high 
temperature limit, c B is zero for n > 6. A model described by resonances of hadron gas in the low 
temperature phase predicts [p(jUs) — p(0)]/T A oc cosh^s/Y) and p„ = \J (n + 2)(n + 1). There- 
fore, both the convergence radiuses in the high temperature and low temperature limits are infinity. 
On the other hand, by using an appropriate scaling ansatz for the free energy at jj, B = 0, one can 
show that C4 will develop a cusp in the 2 flavor chiral limit with rather heavy strange quark mass. 
Hence, cf/cf = p 2 2 should have a peak near the transition temperature and the radius of conver- 
gence may be short near T c when the u, d quark mass is sufficiently small. This implies that the 
distance to the critical values of Hb/T may be estimated by measuring p„ with rather small n. 

The convergence radius p„ have been studied in a 2 flavor simulation with a pion mass of 
m n w 770MeV [pq]. The results of X4/X2 = l^cf /cf are shown by green symbols in FigJ^, where 
Xn = n\c B . The result is consistent with the hadron resonance gas prediction at low temperature 
and with the free gas value at high temperature. However, there is no peak around T c . Since % B 
also increases sharply just below T c , X4/X2 ^ oes not i ncrease near T c although xt itself has a 
pronounced peak. It is interesting to study the behavior of p„ with small u, d quark masses near 
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Figure 3: Quark number density (left) and quark number susceptibility (right) as functions of temperature 
and jlq/T by a simulation with an improved Wilson quark action [^]. Tq is T c at pt q = 0. 



the physical point. The red symbols in Fig § are results of X4/X2 near tne physical mass point in 
2+1 flavor QCD obtained by the RBC-Bielefeld Collaboration [gj. They observed a peak near T c 
beyond the hadron resonance gas value for N t = 4, which may be related to the critical point. 

Moreover, because xl = Xb/T 2 and xf = d{XB/T 2 )/d{}X B /T) at }X B = 0, the figure of X4/X2 
indicates that the baryon number susceptibility (Xb/T 2 ) increases more sharply near the transition 
point as the density is increased when u, d quark masses are small. 

2.4 Study by a Wilson type quark action 

Most lattice QCD studies at finite temperature and density have been performed using stag- 
gered type quark actions with the fourth-root trick of the quark determinant. The theoretical base 
for the fourth-root trick is not confirmed. Moreover, the scaling properties universal to the three- 
dimensional 0(4) spin model, as expected from the effective sigma model, has not been confirmed 
in 2 flavor QCD. Therefore, it is important to carry out simulations adopting different lattice quark 
actions to control and estimate systematic errors due to lattice discretization. 

The WHOT-QCD Collaboration studied finite temperature and density QCD using the clover- 
improved Wilson quark action and the RG-improved Iwasaki gauge action. In contrast to the case 
of staggered quarks, the subtracted chiral condensate shows the scaling behavior with the critical 
exponents and scaling function of the 0(4) spin model for this action [31], and the EoS at \i q = 



have been studied [32]. They performed simulations on 16 3 x 4 lattice along lines of constant 
physics with the mass ratio of pion and rho meson m n /m p = 0.65 and 0.80, and calculated the EoS 
at finite density [Q, [33[ ]. Because the study by a Wilson quark action is more difficult than that by 
staggered quarks in general, some improvements are required. They used a hybrid method of the 
Taylor expansion and the reweighting. Evaluating the quark determinant by the Taylor expansion 
up to 0{pLq), ^ (/J. q ) / ^ (0) was computed. They then obtained the quark number density and the 
susceptibility by numerical differentiations with respect to \i q . The results of the quark number 
density and the susceptibility are plotted in Fig. || as functions of T for each }i q /T . These are quite 
similar to the results obtained by the previous staggered quark simulations. The quark number 
density increases sharply near T c and the slope becomes larger as \l q jT increases. Also, they found 
that a peak in Xq/T 2 appears near T c for large jJ. q /T, suggesting the existence of the critical point. 
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Figure 4: Phase diagram in the (T,ji q ) plane (left) and quark mass dependence of the order of phase transi- 
tions, Columbia plot (right). 




Figure 5: (left) Critical surface in the (m uc i,m s ,jx q ) parameter space predicted by the PNJL model [^]. 
(right) dB4/d((/Xia) 2 ) vs. (/i,a) 2 obtained by a simulation with an imaginary chemical potential [^TJ. L is 
the spatial extent N s . 



3. Critical point at finite density 

In this section, we discuss the critical point terminating the first order phase transition line in 
the (T,ji. q ) phase diagram sketched in Fig. |] (left). The critical point is one of the most interesting 
features that may be discovered in heavy-ion collision experiments. We summarize the current 
discussion on the existence of the critical point in the QCD phase diagram. 

3.1 Quark mass dependence of the critical point 

The order of the phase transition depends on the quark mass for 2+1 flavor QCD. By changing 
the quark mass, the critical point at finite density can be shifted to the low density regime, where we 
can study it by a simulation. The expected nature of the phase transition at pL q = is summarized 
in the right panel of Fig. [| The horizontal axis m u d is the u and d quark masses and the vertical 
axis m s is the strange quark mass. We expect that the phase transition of 2 flavor QCD in the chiral 
limit, (m u d,m s ) = (0,°°), is of second order and that of 3 flavor QCD, (m M ^,m v ) = (0,0), is of first 
order [p4[|. The quenched limit, (m ui j,m s ) = (°°,°°), is also of first order p5| , [36[ ]. The transition 
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for 2 flavor QCD (m s = °°) with finite m uc i is crossover, and 2+1 flavor QCD has a second order 
critical line separating the first order region at small mass and the crossover region at large mass, 
which is shown by the bold red line in Fig. |4] (right). 

We can also discuss the nature of the phase transition at finite density. The left panel of Fig. || 
is a prediction of the critical surface in the (m U( j,m s ,jj, q ) parameter space from the Nambu-Jona- 
Lasinio model with the Polyakov loop (PNJL model) jj37|]. The red lines indicate the critical surface 
which separates the first order and crossover regions. We expect that the first order region becomes 
wider as jj, q increases and the crossover transition at low density changes to be of first order at high 
density for the physical quark masses. 

Mean field argument near the tricritical point 

Let us start with a mean field analysis in the standard sigma model. We discuss the tricritical 
point on the m uc ] = axis in Fig. || (right), at which the second order critical line separates from the 
axis. In the vicinity of the tricritical point at \i q = 0, the effective potential in terms of the chiral 
order parameter a is modeled by the following equation, 

V e n(o) = laa 2 + \ba 4 + ^ca 6 -ho, (3.1) 
Z 4 o 

where we assume c > so that V e ff is bounded from below for large \a\. The coefficients, a,b and 

h may be parameterized as 

? ? T — T c rriE — m s 
a = a t t + a ti jX q , b = b s s + b ti jX q , h = m ud , t = — - — , s= , (3.2) 

where ni£ is m s at the tricritical point. The coefficient b controls the order of phase transition. 
Assuming a symmetry under \i q to — \i q , the leading contribution to b must be pL q at low density. 
Since the effective potential is 0(a 4 ) on the second order critical surface, 

do" 

Solving these equation, we obtain 

The critical surface in the (m U( i,m s ,iJ. q ) space is described by 

c ud m Zi + c s( m E ~ m s) + P^q = 0- ( 3 - 5 ) 
with appropriate constants c u d and c s . The strange quark mass dependence and the \i q dependence 
of the critical light quark mass m c ud around the tricritical point are 

Kd ~ («£ - ™ s ) 5/2 , m c ud ~ (3.6) 

The first equation describes the critical line on the \i q = plane sketched in Fig. f| (right). We 
expect from the first equation that the critical m m i increases very slowly as m s decreases. Similarly, 
the second equation suggests that the chemical potential dependence of the critical surface m^(/i 9 ) 
is also small in the low density region, since the jj. q dependence starts from a term of at = 0. 
The information of the critical surface is important to know the order of phase transition for the real 
world, and the critical surface can be measured near the critical line at \i q = because the study by 
Monte-Carlo simulations is possible in the low density region. 



0, (n = 1,2,3). (3.3) 
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Numerical study of the critical surface 

To investigate the critical surface, some groups performed simulations near the critical quark 
mass at \i q = in QCD with 3 flavors having degenerate quark masses, m u d = m s , and studied the 
\L q dependence of the critical quark mass m c {pL q ). For extrapolating m c {pL q ), an approach on the 
basis of the Taylor expansion in terms of \l q jT j j38| ] and that of the imaginary chemical potential 



[39] have been developed. Moreover, a study of phase-quenched finite density QCD, in which the 



effect from the complex phase of the quark determinant is neglected, has been discussed in []4QfJ . As 
we expect form the mean field argument, the \i q dependence of the critical mass have been found 
to be small in the low density region. 



Recently, an interesting result was obtained by de Forcrand and Philipsen [41]. They studied 
the jU 9 -dependence of the critical quark mass for QCD with 3 flavors of standard staggered quarks 
very precisely and found that the critical line moves towards lighter quark masses as a function 
of pL 2 . They performed simulations with an imaginary chemical potential, \i = \i q a = i^a, on an 
N s x N t = 8 3 x 4 and 12 3 x 4 lattices. Because (detM(ju))* = detM(— /I*) for a complex ju, the 
quark determinant is real if /I is purely imaginary, hence the simulations are possible. In order 
to identify m c (/i,), they computed the fourth order Binder cumulant constructed from the chiral 
condensate, B4 = ((dij/y) 4 ) / ((dxj/y) 2 ) 2 . It has been verified that the critical point belongs to the 
Ising universality class and the value of B4 at T c for m c is the same with that of 3-dimensional Ising 



model, B4 = 1.604, [42]. In the low density region, the analytic continuation from imaginary \x to 



real /I is performed assuming 

B 4 = 1.604 + b l0 (m-m° c ) + b iH 2 + b Q2 H 4 + ■ ■ ■ . (3.7) 

Since bio(m — m®) +Z?oi£i 2 = along the critical line in the leading order, one can estimate the 
curvature of the critical quark mass at /x = by dm c /d{ix 2 ) ~ — boi/bio- Moreover, because B4 
increases when the first order phase transition changes to crossover, bio is positive. Hence, the 
curvature is positive (negative) if b i is negative (positive). The results of ^B4/^((/x ^ •a , ) 2 ) is plotted 
as a function of \L 2 in Fig. || (right), boi is given by bo\ = —Yim^i^dBn/d^jlja) 2 ). From this 
figure, dm c /d(jJ. 2 ) is found to be negative. This result contradicts to the naive expectation. 

In this situation, a simple extrapolation of the critical surface to the physical quark mass point 
is difficult because the first order region in the (m uc i,m s ) plane becomes smaller as \i q increases if 
we do not consider higher order terms of \i q in the Taylor expansion of m c (fi q ). To understand the 
critical surface in the (m u d,m s ,ijL q ) space, studies in a wide range of the chemical potential may be 
necessary. The analytic continuation in the imaginary chemical potential approach is usually based 



on the Taylor expansion of B4 and m c (fl q ), e.g. Eq. (3.7). One of the possible improvements to 



study in the wide range is to use another assumption which is based on a phenomenological model. 



The analytic continuation with various assumptions has been discussed in p3| , 44Q . 

3.2 Reweighting method and Sign problem 

Because the Boltzmann weight is complex at \l q 7^ 0, the Monte-Carlo method is not applicable 
directly. A popular approach to avoid this problem is the reweighting method [45, p6{ ]. We perform 



simulations at [X q = 0, and incorporate the remaining part of the correct Boltzmann weight for finite 
H q in the calculation of expectation values. Expectation values (&) at (jS,/^) are thus computed 
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Figure 6: Histograms of the complex phase 9 [ 12|. Dashed lines are the fit results by Gaussian functions. 



by a simulation at (j6,0) using the following identity: 



'(detAf(/x,)/detAf(0)) 



0,0) 



(detM(ju 9 )/detM(0)) Arf 



(3-8) 



This is the basic formula of the reweighting method. However, because Ae.tM{pL q ) is complex, 
the calculations of the numerator and denominator in Eq. ( |3T8| ) becomes in practice increasingly 
more difficult for larger \i q . If the typical value of the complex phase of the quark determinant 
6 becomes larger than n/2, the real part of e' e (= cos0) changes its sign frequently. Eventually 



both the numerator and denominator of Eq. ( |3.8| ) become smaller than their statistical errors and 



Eq. ( p.8| ) can no longer be evaluated. We call it the "sign problem". 

We define the phase of the quark determinant 6 by the imaginary part of Nf lndetM(^) in the 
framework of the Taylor expansion. In this framework, lndetM(/i 9 ) can be separated into real and 
imaginary parts easily because the even derivatives of lndetM(/i 9 ) are real and the odd derivatives 
are purely imaginary [|24|]. The complex phases 6 are thus given by 



d = ATflm [ln(detM)] = N { £ 



1 



n=0 



d 2,!+1 (lndetM) 2,1+1 
: m— "I' \t) 



{2n+\)\ d{\i q lT) 2n 



(3.9) 



where one must replace A^- in these equations to Nf/4 when one uses a staggered type quark action. 
In Fig. ^, we plot histograms of 6 calculated at the pseudo-critical temperature (/3 = 3.65) for 
jX q /T = 1.0 and 2.0 using the data of the Taylor expansion coefficients up to obtained with 

2 flavors of p4-improved staggered quarks in [|26|]. The application range of the reweighting method 
can be estimated from the histogram. When the width of the distribution is larger than 0(n), the 
phase factor changes its sign frequently. Then, the sign problem happens and the reweighting 
method does not work. Here, it is worth noticing that 6 corresponds to the complex phase of the 
quark determinant however this quantity is not restricted to the range from — % to % because there is 



no reason that the imaginary part of IndetM in Eq. (3.9) must be in the finite range. An interesting 
point which is found from this figure is that these histograms seem to be almost Gaussian functions. 
We fit these data by Gaussian functions, ~ exp(— xd 1 ), with a fit parameter x. The dashed lines in 
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Fig. ^ are the fit results. We find that the histogram of 8 is well-represented by a Gaussian function. 
Such a Gaussian distribution is expected when the system size is sufficiently large in comparison 
to the correlation length due to the central limit theorem. Moreover, the Gaussian distribution of 6 
has been discussed in chiral perturbation theory |^]. 

Once we assume a Gaussian distribution for 6, the problem of complex weights can be avoided 
@. We calculate the expectation value of (Fe' e )^^ =0) , where F = | detM(^ 9 )/detM(0)| iVf 
and G\ detM(/x 9 )/detM(0)| A ' f for the denominator and numerator of Eq. ( |3.8| ), respectively. If 



the operator & is real, the complex phase is given by Eq. ( p.9| ). We introduce the probability 
distribution w as a function of F and 6: 

w(F',d')= J @U8(F' -F)S(6' -6)(detM(0)) Nf e- s z, (3.10) 

where 8 (a) is the delta function. The distribution function itself is defined as an expectation value 
at fl q = 0, however F and 6 are functions of jJ. q /T obtained by a simulation at \i q = 0. 

Since the partition function is real even at non-zero density, the distribution function has the 
symmetry under the change from 6 to —6. Therefore, the distribution function is a function of d 2 , 

e.g., w(d) ~ exp[— (a^Q 2 + a^d A + £?60 6 H )]. If the fluctuations of the phase is small at small \i q , 

the d 2 term gives a dominant contribution. Moreover, as we discussed, the distribution function is 
well-approximated by a Gaussian function: 



w(F,9) « \l ^^w'(F)exp [-a 2 (F)6 2 ] . (3.11) 



We assume this distribution function in terms of 6 when F is fixed. The coefficient ci2{F) is given 
by \/{2a 2 {F')) = (0 2 S(F> -F)) {T ^ =Q) / <S(F' -F)) {T ^_ Q) = (G 2 ) F ,. 
Then, the integration over 6 can be carried out easily, 

( F Wo, =4/" f ^-"'^F = (, W e-(^%^ . (3.12, 

Since 6 is roughly proportional to the size of the quark matrix M, the value of (0 2 ) F becomes 
larger as the volume increases as well as increasing \i q . Therefore, the phase factor in {F{ji. q )e ie ^ 
decreases exponentially as a function of the volume and \i q . However, the operator in Eq. ( |3.12| ) 



is always real and positive for each configuration, hence the expectation value of {F{}X q )e ; ) is 
always larger than its statistical error in this method. Therefore, the sign problem is completely 
avoided if we can assume the Gaussian distribution of 6. 



3.3 Plaquette effective potential 

One of the most primitive approaches to identify the order of the phase transition is to inves- 
tigate the histogram of a typical quantity such as plaquette, Polyakov loop or chiral condensate 
by Monte-Caro simulations. For the case of the plaquette (P), the distribution function, i.e. the 
histogram, is defined by 

W (P') = f @U 8{P'-P) (detM) Nl e 6pN ^ p . (3.13) 



12 



Recent progress in lattice QCD at finite density 



Shinji Ejiri 



For later discussions, we define the average plaquette P as P = —S g / (6j8A^ s i te ), where N s n e = x 
N t . This is a linear combination of Wilson loops for an improved action. If there is a first order 
phase transition point, where two different states coexist at the transition point, the histogram must 
have two peak at two different values of P corresponding to the hot and cold phases. Such studies 
have been done to confirm that the phase transition of SU(3) pure gauge theory is of first order, and 



the double peaked distribution have been obtained at the transition point [35, Bq]. This method is 



equivalent to other methods to identify the order of phase transitions by the Binder cumulant and 
by the Lee- Yang zero [47]. 



The argument of the distribution was extended to the case of finite density QCD by [12]. Here 
and hereafter, we restrict ourselves to discuss only the case when the quark matrix does not depend 
on j6 explicitly for simplicity. The partition function can be rewritten as 

ar&nq) = J R(p,n q ) w (p,p) dp, (3.i4) 

where w(P,/3) is defined in Eq. ( 3.13D at \x q = and R{P,jx q ) is the reweighting factor for finite \i q 
defined by 

R{pu , jmS^-fmm^ (S(P--P)(de»M W/ Me,M(Q)/'.) 
1 S3>Ud(P-P)(AtlM(<S)f< W-P)) m 

This R(P,jJ, q ) is independent of /3 and can be measured at any /3. Here, (•• -}(r ) means the expecta- 
tion value at jj, q = 0. In this method, all simulations are performed at jx q = and the effect of finite 
jjL q is introduced though the operator detM(jx q ) / detM(O) measured on the configurations generated 
by the simulations at jx q = 0. The distribution function for non-zero \i q is R(P,n q )w(P,fi), and thus 
the effective potential is defined by 

Ve ff (P,/3,^) = -ln[R(P,ii q )w(P,P)] =-\aR(P,n q )+V tS {P,p,0). (3.16) 

The shape of the effective potential can then be investigated at \L q ^ once R(P,fi q ) is obtained. A 
schematic illustration of V e s(P) is shown in Fig. [7] (left). 

First, the peak position of the distribution function moves as \i q changes, which is determined 
by solving 

-jp-(.P,P,Hq) = ^-( P '^' ) " ~itf( P >^ = - (3 - 17) 

Then, the effect from \i q to the peak position is the same as that when /3 (temperature) is changed. 
From the definition at jl q = 0, the weight w(P,/$) and the effective potential becomes 

w(P,ftff) = e 6 ^-^ p w(P,p), Veff(P,j3eff,0) = Vdrto/J.O) - 6(/3 eff -/3)AUP, (3.18) 

under a change from j8 to j8 e ff. Hence, the change from j8 to 

jSeft-(^) = + (6N site y l d(lnR)/dP (3.19) 

corresponds (j3,0) to (j3,/i ? ) for the determination of the minimum of V e ff(P). As we will see, 
the slope of In/? is positive. This explains why the phase transition happens when the density is 
increased as well as with increasing temperature (j8). 
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Crossover 



Reweighting factor 
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0.78 0.8 0.82 0.84 0.86 0.88 0.9 0.92 0.94 



Figure 7: (left) Schematic illustration of the effective potential and the reweighting factor, (right) Plaquette 
histogram and the effective potential at \i q = [ |l2[ ] . 
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Figure 8: Reweigiting factor (left) and its curvature (right) as functions of P [[L2|]. 



Moreover, we find from Eq. ( [3.1 8[ ) that the curvature of V e s(P) does not change under the 
change of /3. This means that the curvature of V e ff(P) is independent of j3. The critical value of \i q 
can be estimated by measuring the curvature of the effective potential, 



d 2 V e ff 
dP 2 



( P ,,,) = -^M( fA) + ^( f ,0) = 0. 



(3.20) 



Because the curvature of V r e ff(P, j3,0) at ji q = is positive and the curvature of V e ff(P, j3,/i 9 ) at a 
second order phase transition point is zero, the parameter range where — In R(P,fi q ) has negative 
curvature is required for the existence of the critical point. 



The probability distribution function at non-zero \i q has been calculated in [ ]12| ] using the data 



obtained in [26] with the 2 flavor p4-improved staggered quark action on a 16 3 x 4 lattice. The pion 
mass is m K pa 770MeV, which is heavier than the physical mass. The distribution function w(P) at 
\i q = 0, i.e. the histogram of P, and the effective potential V e ff (P) are given in Fig. [7] (right) for each 
/3. The values of j8 and T /T c are shown above these figures. The potential V e ff(P) is normalized by 
the minimum value for each /3. Because the phase transition is a crossover transition for 2 flavor 
QCD with finite quark mass, the distribution function is always of Gaussian type, i.e. the effective 
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potential is always a quadratic function. The value of the plaquette at the potential minimum 
increases as j8 increases in accordance with the argument of the potential minimum. 

The reweighting factor R(P,jJ, q ) is plotted in Fig. [8] (left). The quark determinant detM(^) is 
estimated using the data of the Taylor expansion coefficients in [p6|]. Higher order terms than jl q 
order are neglected. The application range was selected by checking the truncation error. Because 
the sign problem is serious for the calculation of R(P,/j, q ), the method discussed in Sec. \2_ was 



used to avoid the sign problem. The dashed lines in Fig. |8| (left) are the results that we obtained 
when the effect of the complex phase e' e is omitted. Because the contribution from the complex 
phase is not very large, the error from the approximation to avoid the sign problem may be small. 

To study the existence of a second order phase transition, we discuss the curvature of the 
potential. The result of the curvature, d 2 i\nR)/dP 2 {P,jjL q ), is plotted as solid line in Fig. |8| (right). 
The magnitude of the curvature of \nR becomes larger as jl q /T increases. The dashed line in Fig. || 
(right) is the result of — d 2 (In w)/dP 2 (P) at \l q = 0. This figure indicates that the maximum value of 
d 2 (\nR)/dP 2 (P,n q ) at P = 0.80 becomes larger than -d 2 {\nw)/dP 2 for \i q jT > 2.5. This means 
that the curvature of the effective potential d 2 V e ff/dP 2 vanishes at jJ. q /T ~ 2.5 and a region of P 
where the curvature is negative appears for large pi q /T , corresponding to a double- well potential. 

Further studies are, of course, necessary for the precise determination of the critical point in 
the (T, pi q ) plane, increasing the number of terms in the Taylor expansion of IndetM and decreasing 
the quark mass in the simulation. In particular, the critical value of \L q is sensitive to the quark mass 
as discussed in Sec. 3.1. However, the argument given above suggests the appearance of a first 
order phase transition line at large \l q jT . 



3.4 Canonical approach 

Another interesting approach is to construct the canonical partition function ^c(T,N) by fix- 
ing the total quark number (N) or quark number density (p). Using the canonical partition function, 
we can also discuss the effective potential as a function of the quark number. In this section, we 
denote the grand partition function as 2¥gc{T , pL q ) to distinguish it from the canonical partition 
function 3fc{T,N) explicitly. The relation between 2fQc{T,ji q ) and £?c(T,N) is given by 

3foc{T,ll q ) = f@U (detM(LL q /T)) N <e- s * = £ 3f c (T,N)e N ^ T . (3.21) 
J TV 

Because this equation is a Laplace transformation from i^c to essentially, the canonical parti- 
tion function is obtained from ^Qc{T,jjL q ) by an inverse Laplace transformation. 

In order to investigate the net quark number giving the largest contribution to ^gc(T ', jU 9 ), it is 
worth introducing an effective potential V e ff as a function of N, 

V e n(NJ^ q )^-ln3fc(T,N)-N^ = ^^-N^, Jgc^) =£ e' v ^, (3.22) 

ill N 

where / is the Helmholtz free energy. Using the effective potential for the quark number, the 
argument of the nature of phase transition is possible as well as the effective potential for the 
plaquette, and the physical meaning of this effective potential is clearer than that of the plaquette. 
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Figure 9: Phase structure in the (T,p) plane and the behavior of fl*/T as a function of p. 



If there is a first order phase transition region, we expect that this effective potential has minima 
at more than one value of N. At the minima, the derivative of V e ff satisfies 

^CW.r.M,)— *^(T,*0-*-0. (3.23) 

Hence, in the first order transition region of T, we expect d (In 3fc) I dN (T ,N) = —flq/T takes the 
same value at different N. Here, H*(T,N) is the chemical potential which gives a minimum of the 
effective potential at (T,N) and becomes \i q in the thermodynamic limit because the potential is 
minimized in the large volume limit. 

The phase structure in the (T,p) plane and the expected behavior of jU*/r are sketched in 
the left and right panels of Fig. [| respectively. The thick lines in the left figure are the phase 
transition line. We expect that the transition is crossover at low density and becomes of first order 
at high density. Since two states coexist on the first order transition line, the phase transition line 
splits into two lines in the high density region, and the two states are mixed in the region between 
two lines. The expected behavior of fl* along the lines A and B are shown in the right figure. 
When the temperature is higher than the temperature at the critical point T pc (line A), jj.* increases 
monotonically as the density increases. However, for the case below T cp (line B), this line crosses 
the mixed state. Because the two states of pi and P2 are realized at the same time, fi* does not 
increase in this region between pi and P2. 



Glasgow method [ }45| , |48Q has been a well-known method to compute the canonical parti- 
tion function. A few years ago, such a behavior at a first order phase transition was observed by 
Kratochvila and de Forcrand in 4 flavor QCD with staggered fermions [ |Tl| ] calculating the quark 
determinant by the Glasgow algorithm on a small lattice. Also, Alexandru et al. [|^] proposed 
a method to perform simulations with canonical ensemble directly, and recently the method was 



tested for 2 and 4 flavor QCD [14]. Moreover, a method based on a saddle point approximation 
was proposed by [13]. By this approximation, the computational cost is drastically reduced and the 
first order like behavior was observed for 2 flavor QCD. We explain these recent developments. 

Inverse Laplace transformation 



From Eq. (3.21), the canonical partition function can be obtained by an inverse Laplace trans- 
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Nf=4, Wilson gauge and fermion action, 6 3 x4 




Nf=2, Wilson gauge and fermion action, 6 3 x4 
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Figure 10: Baryon chemical potential as a function of baryon number for Nf = 4 (left) Nf = 2 (right) 
obtained by simulations with canonical ensemble [114], 



formation [50, 51], 



iT c (r,A0 = |- f /3 e-^/ r +*/ r )jGc(r,Mo + iMO , (3-24) 



where /lo is an appropriate real constant and /x,- is a real variable. Note that J?oc(r, ju 9 + 2%iT /3) = 
3?Gc(T,n q ) [p2|]. The grand canonical partition function can be evaluated by the calculation of the 
following expectation value at ju 9 = 0. 



3bc(r,jUg) = l 

^Gc(r,0) ^gc 



V detM(0) 



detM(0) Nt e 



V detM(0) 



AT f \ 



(7-^=0) 
(3.25) 



If one adopts /^o = 0, the chemical potential in Eq. ( |3.25[ ) is purely imaginary and the quark 
determinant is real. Therefore, the calculation of the quark determinant at i^u, is easier than the 
calculation at real \L q . Performing a simulation at \i q = 0, the standard staggered quark determinant 
can be computed for any chemical potential at the cost of diagonalizing a matrix modified from the 
staggered fermion matrix by using the Glasgow algorithm [48]. Because Eq. Q3.24 ) is actually a 
Fourier transformation, if one calculates the Fourier coefficients of detM(/jU,-/r) for each configu- 
ration, the canonical partition function 3fc(T,N)/ Sqc{T,0) can be computed by averaging them 
over the configurations. 

Kratochvila and de Forcrand calculated —d (In Sc.) I dN(T,N) as a function of the baryon num- 
ber Nb = N/3 with a fixed temperature (j8 ) performing a simulation of 4 flavor QCD with staggered 



quarks on a 6 3 x 4 lattice [11]. The phase transition is of first order even at \i q = for their simu- 
lation of 4 flavor QCD. Their result of — d (In Sq) / dN(T ,N) shows an S-shape in the temperature 
range below T c , as we discussed above for a first order phase transition. 

Simulations with canonical ensemble 



By replacing the order of the integrals of the inverse Laplace transformation Eq. ( p.24[ ) and the 
path integral in Sqc> the canonical partition function is written by 



S C (T,N) = J @Ue- s sdet N M Nf , 



(3.26) 
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Figure 11: Chemical potential vs. quark number density for Nf = 2 with a saddle point approximation [pj|]. 
where det^M^ is defined by 



det/vM' 



Nf 



1 

2n 



2k 



-iNm/T 



&etM{iHi/T) Nt d 



(3.27) 



The method to perform Monte-Carlo simulations with this partition function has been proposed 
by [||, ||]. Although det^M^ is complex, we obtain det N M Ni = (det^ N M Nl )* and 3T C {T,N) = 
3?c(T, —N) from a symmetry under the replacement from \i q to — \i q . Using these properties, we 
can rewrite the partition function as 



(3.28) 



Then, the Boltzmann weight is real and the Monte-Carlo method is applicable. One can generate 
configurations according to the weight exp(— S g )Re(detNM Nf ), and det^M^ is computed numeri- 
cally with an approximation, 



det/vM^— £ e- 2 ^ w /^detM(27r/;/^V dis ) A ' f . 

A'dis y = o 



(3.29) 



This approximation is applicable for A^i s S> A^. 

The ^QCD collaboration (Kentucky group) performed simulations with the canonical ensem- 
ble on a 6 x 4 lattice for 2 and 4 flavor QCD with Wilson quarks using this method and the 



preliminary results are presented in this conference [14]. The results of the chemical potential 
Mb I t = ~d (In Jc) /d (N/3) are shown in Fig. [K)| for 4 flavor QCD (left) and 2 flavor QCD (right). 
Their result of 4 flavor QCD shows an S-shape in the temperature range below T c , suggesting a first 



order phase transition. This result is consistent with the previous result Ql 1Q. On the other hand, the 
chemical potential increases monotonically for 2 flavor QCD. This result suggests that the phase 
transition is crossover at the temperature they investigated. 

Saddle point approximation 

However, the studies by above-mentioned two methods need much computational cost and are 
difficult except on a small lattice with present day computer resources. To reduce the computational 
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cost, a method based on a saddle point approximation has been proposed in [13]. If one selects a 



saddle point as jito in Eq. ( p.24| ). The information which is needed for the integral is only around 
the saddle point when the volume is sufficiently large. Moreover, if we restrict ourselves to study 
the low density region, the value of detM(n q /T) near the saddle point can be estimated by the 
Taylor expansion around jx q = 0. The calculations by the Taylor expansion are much cheaper than 
the exact calculations and the studies using large lattices are possible. Also, the truncation error 
can be systematically controlled by increasing the number of the expansion coefficients. 



We perform the integral in Eq. (3.24) by a saddle point approximation. We denote the quark 
number density in a lattice unit and physical unit as p = N/N 3 and p/T 3 = pN 3 , respectively. 
We assume that a saddle point zo exists in the complex \l q /T plane for each configuration, which 
satisfies D'(z.o)-p = 0, where (detM(z)/detM(0)) A ' f = exp[^ 3 D(z)] and D'(z) = dD(z)/dz, We 
then perform a Taylor expansion around the saddle point and obtain the canonical partition function, 



&c(T,pV) = —%> GC (T,0)n ^xp 







/ , eXP 


H 


J-k/3 





pzo 



iTGc(r,0)(exp[V(D(zo)-pzo)]e 



-ia/2 _ 



V\D"(zo)\. 



dx 



(3.30) 



(7-^=0) 

Here, D"{z) = d 2 D(z)/dz 2 , V = N 3 and D"{z) = \D"(z)\e ia . We chose a path which passes the 
saddle point. Higher order terms in the expansion of D(z) become negligible when the volume V 
is sufficiently large. 

We calculate the derivative of the effective potential with respect to N or p. Within the frame- 
work of the saddle point approximation, this quantity can be evaluated by 



T 



1 d In %c(T,pV) 



zo ew[V(D( Z0 )-pzo)]e- ia / 2 



V\D"(zo)\ 



V 



dp 



exp[V (D(z )-pZo)]e- ia / 2 



l 



(3.31) 



This equation is similar to the formula of the reweighting method for finite \i q . The operator in the 
denominator corresponds to a reweighting factor, and jU*/r is an expectation value of the saddle 
point calculated with this modification factor. 

The derivative of In Jfc was computed in [|13|] using the data obtained in [ 26] with the 2 flavor 
p4-improved staggered quark action, m n rj 770MeV. Because the modification factor is a com- 
plex number, this calculation suffers from the sign problem. To eliminate the sign problem, the 



approximation discussed in Sec. 13^ was used. If one assumes that the distribution of the complex 
phase is well-approximated by a Gaussian function, the complex phase factor e' 6 can be replaced 
by exp[— (8 )/2]. Moreover, the quark determinant was estimated by the Taylor expansion up 



to 0{ji. q ). Because the calculation of Eq. ( 3.31 ) is similar to the calculation by the reweighting 



method, the configurations which give important contribution are changed by the modification fac- 



tor. To avoid this problem, the multi-/3 reweighting method [54] was used. By this method, the 
important configurations are automatically selected among all configurations generated at many 
simulation points of /3 . 



The result of \l* q jT is shown in Fig. 11 as a function of p/T for each temperature T/T C (J3). 
The dot-dashed line is the value of the free quark-gluon gas in the continuum theory, p/T 3 = 
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Ni[{jl q /T) + (l/7T 2 )( i u 9 /r) 3 ]. From this figure, we find that a qualitative feature of H q /T changes 
around T /T c ~ 0.8, i.e. ;U*/T increases monotonically as p increases above 0.8, whereas it shows 
an S-shape below 0.8. This means that there is more than one value of p/T 3 for one value of 
/i*/r below T/T c ~ 0.8. This is a signature of a first order phase transition. Although some 
approximations are used, the critical value of jJ.*/T is about 2.4, which is roughly consistent with 
the critical point estimated in Sec. by calculating the effective potential of the plaquette using 
the same configurations, (T /T c ,^i q /T) « (0.76,2.5). The difference between these two results may 
be a systematic error. Further studies are necessary to predict the critical point quantitatively, but 
we find that the canonical approach is useful to study the phase structure at finite density. 



4. Summary 

We reviewed recent studies of lattice QCD at finite density. Remarkable progress in the study 
of the equation of state was obtained. The MILC Collaboration and the RBC-Bielefeld Collabo- 
ration performed simulations with quark masses near the physical point using improved staggered 
quark actions and calculated thermodynamic quantities along lines of constant entropy per baryon 
number in the low density region using the Taylor expansion method. Because experimental re- 
sults obtained in heavy-ion collisions can be well-explained by a perfect fluid model, the isentropic 
equation of state is needed for the analysis of the experimental data. Moreover, fluctuations of 
hadron numbers at finite density are important, which can be measured in the experiments. If there 
is a critical point, the fluctuation of baryon should be large around that point. Recent simulations 
show that the baryon number fluctuation makes a peak near T c at finite \i q and the fluctuation be- 
comes larger when the quark masses are decreased to the physical point. Confirmations by other 
quark actions are also important. Simulations with a Wilson type quark action were performed by 
the WHOT-QCD Collaboration and the large fluctuation at finite \i q was discussed. 

The chemical potential dependence of the critical line in 2+1 flavor QCD with quark masses 
(m u d,m s ) was studied at low density to understand the grovel structure of the critical surface in the 
{m u d,m s ,pL q ) parameter space. The \L q dependence of the critical quark mass is found to be small in 
the low density regime and the current result of dm c / d (/j.^) is slightly negative at m uc i = m s . Further 
studies in a wide range of the parameter space are important to understand the phase structure. 

Some methods to investigate finite density QCD beyond the low density region were also 
discussed. A method based on the investigation of an effective potential as a function of the average 
plaquette was proposed introducing an approximation to avoid the sign problem, and the existence 
of the critical point at finite density was suggested by a simulation with improved staggered quarks. 
Moreover, it was found that interesting information about the order of phase transitions at finite 
density is obtained by constructing the canonical partition function for each quark number. 
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